clear; clc;

%% 参数固定
a = 1;
c = 11.5;

% 初始条件
x0 = 0.4; y0 = 0.8; z0 = 0.6;

% 李雅普诺夫指数计算参数
N_LE = 1000;
transient_LE = 200;
delta = 1e-8;

% 分岔图计算参数
N_bif = 1200;
transient_bif = 200;
max_pts_per_b = 100;

%% b取值范围
b_list = linspace(2, 8, 500);

LEs_b = zeros(length(b_list), 3);
X_cell = cell(length(b_list), 1);  % 用 cell 存储分岔数据，避免 parfor 赋值错误

parfor i = 1:length(b_list)
    b = b_list(i);
    
    % 计算李雅普诺夫指数
    LEs_b(i,:) = MLE(a, b, c, x0, y0, z0, N_LE, transient_LE, delta);
    
    % 分岔图计算
    x = zeros(1, N_bif);
    y = zeros(1, N_bif);
    z = zeros(1, N_bif);
    x(1) = x0; y(1) = y0; z(1) = z0;
    
    for k = 1:N_bif-1
        x(k+1) = a * sin(b * z(k)) * sin(c / max(abs(x(k)),1e-6));
        y(k+1) = a * sin(b * x(k+1)) * sin(c / max(abs(y(k)),1e-6));
        z(k+1) = a * sin(b * y(k+1)) * sin(c / max(abs(z(k)),1e-6));
    end
    
    x_vals = x(transient_bif+1:end);
    sample_indices = round(linspace(1, length(x_vals), max_pts_per_b));
    x_sampled = x_vals(sample_indices);
    
    X_cell{i} = [repmat(b, max_pts_per_b, 1), x_sampled'];
end

% 循环外合并分岔数据
X_plot_b = vertcat(X_cell{:});

%% 绘图
figure;

% 李雅普诺夫指数图 (Fig.1d)
subplot(1,2,1);
plot(b_list, LEs_b(:,1), 'r-', 'LineWidth', 1.5); hold on;
plot(b_list, LEs_b(:,2), 'b-', 'LineWidth', 1.5);
plot(b_list, LEs_b(:,3), 'g-', 'LineWidth', 1.5);
xlabel('b');
ylabel('LEs');
xlim([1 8]);
title('Fig.1(d) 李指数图');
legend('LE1','LE2','LE3','Location','Best');
grid on;

% 分岔图 (Fig.1e)
subplot(1,2,2);
plot(X_plot_b(:,1), X_plot_b(:,2), '.', 'MarkerSize', 1);
xlabel('b');
ylabel('x');
xlim([1 8]);
title('Fig.1(e) 分岔图');
grid on;
